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Analysis and denoising of Cosmic Microwave Background (CMB) maps are per- 
formed using wavelet multiresolution techniques. The method is tested on 12°. 8 x 12°. 8 
maps with resolution resembling the experimental one expected for future high resolu- 
tion space observations. Semianalytic formulae of the variance of wavelet coefficients 
are given for the Haar and Mexican Hat wavelet bases. Results are presented for the 
standard Cold Dark Matter (CDM) model. Denoising of simulated maps is carried out 
by removal of wavelet coefficients dominated by instrumental noise. CMB maps with 
a signal-to-noise, S/N ~ 1, are denoiscd with an error improvement factor between 3 
and 5. Moreover we have also tested how well the CMB temperature power spectrum 
is recovered after denoising. We are able to reconstruct the Ct's up to I ~ 1500 with 
errors always below 20% in cases with S/N ^ 1. 
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1 INTRODUCTION 

Future CMB space experiments will provide very detailed 
all-sky maps of CMB temperature anisotropics; NASA MAP 
Mission (Bennett et al. 1996) and the ESA Planck Mission 
(Mandolesi et al. 1998; Puget et al. 1998). The high sen- 
sitivity of these experiments will result in unique data to 
constrain fundamental cosmological parameters. Moreover, 
future CMB maps will allow to distinguish between compet- 
ing theories of structure formation in the early universe and 
will provide very fruitful data on astrophysical foregrounds. 

The cosmological signal in CMB maps is hampered by 
instrumental noise and by foreground emissions. Therefore, 
a necessary step in analysing CMB maps is to separate the 
foreground emissions from the CMB signal. Several linear 
and non-linear methods have already been tested on simu- 
lated data (Bouchet, Gispert & Puget 1996; Tegmark & Efs- 
tathiou 1996; Hobson et al. 1998a,b). An alternative method 
can be one based on wavelets. Wavelets are known to be 
very efficient in dealing with problems of data compression 
and denoising. Development of wavelet techniques applied 
to signal processing has been very fast in the last ten years 
(see Jawerth & Sweldens 1994 for an overview). These tech- 



niques have already been applied to a variety of astrophysi- 
cal problems. For example, regarding cosmology, Slezak, de 
Lapparent & Bijaoui (1993) have applied wavelet analysis to 
the detection of structures in the CfA redshift survey. They 
have also been introduced to study the Gaussian character 
of CMB maps (Pando et al. 1998, Hobson et al. 1998). A 
study using spherical Haar wavelets to denoise CMB maps 
has just appeared (Tenorio et al. 1999). 

We consider small patches of the sky where a flat 2-D 
approach is valid. We apply wavelet multiresolution tech- 
niques, known to be computationally very fast taking only 
O(N) operations to reconstruct an image of N pixels. In 
the 2-D flat wavelet analysis a single scale and two transla- 
tions are usually introduced, where the basis is generated by 
4 tensor products of wavelets and scaling functions. There- 
fore, three detail images plus an approximation image ap- 
pear at each level of resolution. Different wavelet bases are 
characterized by their location in space. The bases consid- 
ered in this work are Mexican Hat, Haar and Daubechies. 
The first one is the most localized though, as opposed to the 
other two, it does not have a compact support. As a first ap- 
proach to the application of these techniques to CMB data 
we only consider maps with cosmological signal plus instru- 
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Figure 1. Variance of diagonal (solid lines) and horizon- 
tal/vertical (dotted lines) detail wavelet coefficients Cd,Ch versus 
scale R. A standard CDM cosmological model is assumed. Thin 
lines outline the result obtained for the Haar wavelet system and 
the thick lines correspond to the Mexican Hat wavelet basis. 



mental Gaussian noise. CMB experiments are contaminated 
by galactic (dust, free-free and synchrotron emission) and 
extragalactic foregrounds (infrared and radio galaxies, S-Z 
effects from clusters,...) in addition to noise. Obviously in 
view of this, the contents of this paper only cover a small 
fraction of the work to be done. Our aim is to shed light 
on the wavelet characterization of the different components, 
the late phase being to build up a wavelet-based framework 
to disentagle all of them. In this line, a Bayessian method 
(maybe incorporating entropy or other constraints) dealing 
with wavelet components at different scales and integrating 
the different channels will be the final goal. Knowing the 
efficiency of wavelets in removing noise as shown by many 
works in other fields, the first step to take in the application 
of wavelet techniques to the CMB is to study the denoising 
of temperature maps. 

The outline of the paper is as follows. A theoretical con- 
tinuous wavelet analyis of CMB data is presented in Section 
2. General semianalytical formulae are given for the variance 
of the detail wavelet coefficients as a function of the temper- 
ature power spectrum Ct . Section 3 introduces the discrete 
wavelet technique that is applied to denoising of simulated 
CMB maps in Section 4. Discussion and conclusions are pre- 
sented in Section 5. 



2 CONTINUOUS WAVELET ANALYSIS 

2.1 One-dimensional Transform 

The Fourier transform is a powerful tool in many areas 
but in dealing with local behaviour shows a tremendous 
inefficiency. For instance, a large number of complex ex- 
ponentials must be combined in order to produce a spike. 
The wavelet transform solves this problem, introducing a 
good space-frequency localization. It is conceptually sim- 



ple and it constitutes a fast algorithm. Let if)(x) be a one- 
dimensional function satisfying the following conditions: a) 
f_ dxip(x) = 0, b) r* 5 dxip 2 {x) = I and c) Cp = 

J_ dk\k\~ 1 ip(k) < oo, where tf)(k) is the Fourier trans- 
form of tp(x). So, according to condition a), the wavelet must 
have oscillations. Condition b) is a normalization and c) rep- 
resents an admissibility condition in order to reconstruct a 
function f(x) with the basis if) (see equation (2) for such a 
synthesis) . 

We define the analyzing wavelet as *$>(x;R, b) = 
R~ ' "i/H^Tp), dependent on two parameters: dilation (_R) 
and translation (6). It operates as a mathematical micro- 
scope of magnification _R _1 at the space point b. The wavelet 
coefficients associated to a one-dimensional function f(x) 
are: 



w(R,b) = / dxf(x)V(x;R,b) 



(1) 



It is clear from the above definition that such coefficients 
represent the analyzing wavelet at x for a delta distribu- 
tion peaked at this point, i.e. for f(x) = S(x — x ). For 
R — 1, w(R, b) is the convolution of the function / with the 
analyzing wavelet if). 

The reconstruction of the function / can be achieved in the 
form 



/(Z) = (2t:C^- 1 



dRdb R~ 2 w(R, b) *(z; R, b) 



(2) 



Examples of wavelet functions are: i) Haar, if> — 1( — 1) for 



< x < 1/2 (1/2 < x < 1), ii) Mexican hat, if) ■■ 



(9tt)V 



174(1- 



2 /2 



2.2 Two-dimensional Transform 

Regarding the two-dimensional case, we introduce a one- 
dimensional scaling function <f> normalized in the form: 
\_ dx(j>(x) — 1. Examples of scaling functions are: i) 
Haar, <f> = 1 (0) for < x < 1 (x < 0, x > 1), ii) 



Mexican Hat, 



(9 7r ) 1 /4 e 



-x^/2 



The analyzing scaling 



&(x;R,b) = R i// 0( 2 ^), allows to define details of an im 



age, f(x), with respect to the tensor products 

T d (x; R, b) = ®(xi;R, fei)*(x 2 ; R, b 2 ) , (3) 

F h (x;R,b)=®(xi;R,bi)*(x2;R,b 2 ) , (4) 

r v (x ■ R, b) = 9(xy,R, bi)$(xr,R, h) . (5) 

The diagonal, horizontal and vertical wavelet coefficients are 
defined by (a = d, h, v) 



w a (R,b) 



dx f(x) T a (x; R,b) 



(6) 



Scaling functions act as low-pass filters whereas wavelet 
functions single out one scale. Therefore, detail coefficients 
provide local information about symmetrical (diagonal) and 
elongated/filamentary structure (vertical and horizontal). 

Let us now assume an homogeneous an isotropic ran- 
dom field f{x), i.e. the correlation function C(r) =< 
f(x)f(x + r) >, r = \r\, where <> denotes an average 
value over realizations of the field. The Fourier transform 
of the field /(ft) satisfies < f{k)f{k') >= P(k)5 2 {k - ft'), 
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Figure 2. rms deviation of wavelet detail coefficients obtained from CMB maps (signal maps, dashed-dotted lines), CMB plus noise 
maps (signal plus noise maps, solid lines) and pure noise (dashed lines) are presented in left panels. Right panels show the ratio of the 
rms deviation of the detail coefficients from signal maps divided by the rms deviation of the detail coefficients from noise maps. Top 
panels correspond to simulated maps with S/N = 0.7; bottom panels correspond to S/N = 2. 



where k = \k\ and P(k) is the power spectrum (the Fourier 
transform of C(r)). In this case we can calculate the corre- 
lation and variance of the wavelet coefficients: C a {r; R) =< 
w a (R,b)w a {R,b + r) >,<Ja(R) = C a {0; R) and we find the 
following equations 



C(0) = a 



C?l dRR~ 3 a 2 a (R) 



(2tt) 2 / dkk- 2 \Y a \ 2 {k) 



(7) 



where f a (k) is the Fourier transform of RY a . 

On the other hand, we calculate the Fourier transform 
of the wavelet coefficients w a (R, b) with respect to the b 
parameters: 

< w a (R,k)w a >(R',k') >=w aa ,(R,R';k)S 2 (k-k') , (8) 

WW = (2-K) 2 RR'P{k)F* a {Rk)r a ,{R'k) , (9) 

that allows us to get the detail wavelet variances as 

(10) 



a 2 a {R) = / dkP(kR~ 1 )\f a (k)\ 
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The diagonal variance corresponds to the tensor prod- 
uct of two one-dimensional wavelets. If ^(fc)! 2 is a function 
strongly peaked near k ~ 1 then a 2 (R) ~ P(k ~ -R -1 ), 
taking into account the normalization of the wavelet func- 
tion, that allows an estimation of the power spectrum in 
terms of the diagonal component. This is what happens 
for the Mexican hat: |i/>(fc)| 2 oc k 4 e~ k , with a maximum 
at k = 2" 1 ' 2 , whereas the Haar wavelet is not localized 
in Fourier space: |i/i(fc)j 2 oc (k / A)~ 2 sin A (k / A) . We can also 
deduce that Cu = C v and a 2 — a 2 taking into account 
the symmetry of the equations. Moreover, the tempera- 
ture power spectrum P(k) can be obtained from the detail 
wavelet power spectrum w aa (R,R;k) as follows 



/'(/-■) = — j -^ / ,l()u:„An.lhl,;i) 



1 

cv~ 



n — (cos 9 , sin 9) . (11) 

For the Haar and Mexican wavelets we can calculate: 



HAAR: \r d \ 2 



1 iife " 
(2tv) 2( 4 ' 



r . fcl . fe2,4 

[ Sm _ Sm _] 



|r„| = 



(27T) 



1 f hb , 

2 V a - 



r - k 2 ,4,l . 2, 
sin — — sm ki 
4 4 



(12) 
(13) 



MEXHAT: \T d \~ = —(kika)'e 

97T 



|f,| 



4 L. 2 -k 2 

-z-k 2 e 
■iix 



(14) 



where k 2 = k\ + fcf ana f« can be obtained from f h, swap- 
ping fcj and k%- The variance of the detail wavelet coeffi- 
cients for the Haar and Mexican Hat systems, assuming the 
standard CDM model, is presented in Figure 1. As one can 
see the acoustic peaks can be clearly noticed, being more 
pronounced for the Mexican Hat basis. This last result is a 
consequence of being a more localized wavelet system. For 
a more detailed discussion see Sanz et al. 1998, 1999. 



3 DISCRETE WAVELET ANALYSIS 

3.1 One-dimensional Multiresolution Analysis 

An orthonormal basis of L 2 (SR) can be constructed from a 
wavelet ip through dyadic dilations j and translations k 



ip 3 , k {x) = 2 j/2 4>(2 j x-k) 



(15) 



In addition, a scaling function <j> can be defined associated 
to the mother wavelet ip. Such a function gives rise to the 
so called multiresolution analysis. A multiresolution analysis 
of L 2 (5R) is defined as a sequence of closed subspaces Vj of 
L 2 (5R), j € Z. Properties can be seen in Ogden (1997). 

Subspaces Vj are generated by dyadic dilations and 
translations of the scaling function <p (this function forms 
an orthonormal basis of V , {<f> ,k{x) = 4>(. x ~ &)})• More- 
over each Vj can be expressed as the orthogonal sum Vj = 
Vj-i (B Wj-i, where Wj-i is created from wavelets ipj^i^- 
Taking into account the properties of the scaling function, 
together with this last expression, we can construct approxi- 
mations at increasing levels of resolution. These approxima- 
tions are linear combinations of dilations and translations of 
a scaling function <f>. The difference between two consecutive 
approximations, i.e. the detail at the corresponding resolu- 
tion level, is given by a linear combination of dilations and 
translations of a wavelet function ip. 



then (r a (af; j, k)T a ,(x;j' ,k')) - S aa >8jjiS%p, where () de- 
notes the scalar product in L 2 (K 2 ). If we define the discrete 
wavelet coefficients associated to any detail by the equation 
(6) 

w a (j,k) = / dxf(x)r a (x;j,k) , (18) 

we can thus reconstruct the image with all the details 



/(*) = 






_w a (j,k)T a (x;j,k) 



aj,* 



(19) 



In particular, we get the following expression for the second- 
order moment of the image 

(/ 2 0?))=V w 2 a {j,k), (20) 

' 'a,j,k 

that expresses how the energy of the field is distributed lo- 
cally at any scale and detail. 

For a finite image, Rmax X Umax, in order to re- 
construct it we must add to equation (19) an approx- 
imation w a (k)r a (x;k) with F a (f;fc) = 3>(xy, Rmax, ki) 
$(x2;R m ax,k 2 ) and w a (k) = J dx f(x) Y a (x, k), represent- 
ing the field at the lower resolution. If f(x) represents the 
temperature fluctuation field then the variance is given by 
< (AT/T) 2 >= ((AT/T) 2 )/N P , being N p the number of 
pixels. 

The orthonormal basis that we are going to use are the 
standard Daubechies TV (Haar corresponds to TV = 1), that 
has been extensively used in the literature because of their 
special properties: they are defined in a compact support, 
have increasing regularity with TV and vanishing moments 
up to order TV — 1 (Daubechies 1988). On the contrary, the 
Mexican Hat wavelet is not defined in a compact support 
and it is not appropriate for this multiresolution analysis. 

For discrete wavelet analysis of the CMB maps we use 
the Matlab Wavelet Toolbox (Misiti et al. 1996). This tool- 
box is an extensive collection of programs for analyzing, de- 
noising and compressing signals and 2-D images. Discrete 
Wavelet decomposition is performed as described above to 
obtain the approximation and detail coefficients of the 2-D 
CMB maps at several levels. 



3.2 Two-dimensional Multiresolution Analysis 

The analysis performed in this work assumes equal dila- 
tions in the 2 dimensions involved. At a fixed level of reso- 
lution, subspaces in a 2-D multiresolution analysis are the 
tensor products of the corresponding one-dimensional ones 
Vj+i = Vj+i <8> Vj+i.The 2-D basis is therefore built by the 
product of two scaling functions (approximation), the prod- 
uct of wavelet and scaling functions (horizontal and vertical 
details) and the product of two wavelets (diagonal details): 

Vj+i = (Vj®Wj)®(Vj®Wj) (16) 

= (Vi®V5)et(V5®w r i )e(w>®v i )©(w i ®Wj)] .(17) 

Horizontal, vertical and diagonal detail coefficients represent 
the variations in these directions relative to a weighted av- 
erage at a lower resolution level (given by the approximation 
coefficients) . 

A discrete orthonormal basis, T a (x;j,k), can be de- 
fined by setting R — 2~ 3 and b — 2~ 3 k in equations (3-5), 



4 DENOISING OF CMB MAPS 

Future CMB space experiments will provide maps with res- 
olution scales of few arcminutes. In this work we analyze 
simulated maps of 12.8 x 12.8 square degrees with pixel size 
of 1.5 arcmin. Simulations are made assuming the standard 
CDM, n = 1 and H = 50 km/sec/Mpc. The maps are fil- 
tered with a 4'. 5 FWHM Gaussian beam to approximately 
reproduce the filtering scale of the High Frequency Chan- 
nels of the Planck Mission. Simulated maps have a rms 
signal of AT/T — 3.7 x 10~ 5 . Gaussian noise is added to 
these maps at different S/N levels between 0.7 and 3. A 
non-uniform noise is also considered to account for the non- 
uniform sampling introduced in satellite observations. As an 
extreme case we have assigned the signal-to-noise at each 
pixel from a truncated (at the 2<r level) Gaussian distribu- 
tion with a mean value of 2 and a dispersion of 0.5. We 
use the set of Matlab Wavelet 2-D programs with the cor- 
responding graphical interface to analyze and denoise those 
maps. Suitable bases of wavelets are studied. Daubechies 4 
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Figure 3. Mean value (solid line) and la error (dashed-dotted lines) of the absolute value of the relative errors, ACf/Cf. Top-left panel 
corresponds to S/N = 1.0, top-right to S/N = 2, bottom-left to S/N = 2 with non-uniform noise and bottom-right to S/N = 3. 



wavelets are the ones used in this analysis. No significant 
changes are observed when the analysis is carried out using 
other higher order Daubechies bases. On the other hand, 
the Haar system is not so efficient for denoising CMB maps 
since it produces reconstruction errors much larger than us- 
ing high order Daubechies systems. 

First of all, three wavelet decompositions are performed 
obtaining wavelet coefficients corresponding to the CMB 
original map, to the signal plus noise map and to the pure 
noise map. Decompositions are carried out up to the fourth 
resolution level. Denoising of the signal plus noise maps is 
based on subtraction of certain sets of coefficients affected 
by noise. White noise is the most common in CMB exper- 
iments. The dispersion of wavelet coefficients of that type 
of noise is constant as can be seen from equation 10. On 
the contrary CMB detail wavelet dispersions go to zero as 



R goes to zero. Therefore first level wavelet coefficients are 
dominated by noise and then, for a given signal plus noise 
map, it is possible to know the noise and consequently the 
CMB wavelet coefficient dispersions at all levels. CMB maps 
produced by typical experiments with a ratio between an- 
tenna and pixel size of « 3 will have wavelet coefficients 
containing the relevant information on the signal at level 3 
and above. As shown below, level 3 is the critical one to per- 
form denoising as the noise can still be at a level comparable 
to the signal. Figure 2 shows rms deviations and correspond- 
ing ratios for two simulations with S/N = 0.7 and S/N — 2. 
Detail coefficient numbering corresponds to the three direc- 
tions diagonal, vertical and horizontal at the three consec- 
utive levels, i.e., numbers 1,2,3 correspond to diagonal, ver- 
tical and horizontal cofficients at the first resolution level, 
4,5,6 to the second level coefficients in the same order and 
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Figure 4. Absolute value of the relative errors, ACi/Ce, of the CMB power spectrum obtained from signal-plus-noisc maps (solid lines), 
wavelet denoised maps (short dashed lines) and Wiener denoised maps (dashed lines). Top-left panel corresponds to S/N = 1 (wavelet 
denoised maps removing all coefficients at levels 1, 2, 3d and 3h is included as long dashed lines), top-right to S/N = 2, bottom-left to 
S/N = 2 with non-uniform noise and bottom-right to S/N = 3. 



7,8,9,10,11,12 to levels 3 and 4 respectively. As it can be 
seen, the first two levels are entirely dominated by noise as 
pointed out before. Therefore, all these coefficients can be 
removed to reconstruct a denoised map. This is equivalent 
to using a hard thresholding assuming a threshold above all 
these coefficients. On the other side, level 4 is completely 
dominated by the CMB signal and is left untouched. Ratios 
between rms deviations of the signal and noise maps at the 
third resolution level are not always clearly dominated by 
noise or signal. Ratios of w 1 are treated with a soft thresh- 
olding technique (in practice we consider ratios in the range 
0.3 — 1.5 though changes in this interval do not significantly 
affect results). Soft thresholding consists of removing all co- 
efficients with absolute values smaller than the threshold de- 
fined in terms of the noise dispersion (cr n ). Coefficients with 



absolute values above the defined threshold are rescaled by 
subtracting the threshold to the positive ones and adding it 
to the negative ones. To define these thresholds we use the so 
called SURE thresholding technique introduced by Donoho 
& Johnstone 1995. This technique is based on finding an es- 
timator of the signal that will minimize the expected loss or 
risk defined as the mean value of (1/N P ) ^2-- i(^d< ~ Ti) 2 , 
where Ti is the temperature at pixel i in the original signal 
map and Td i is the estimator at pixel i (temperature in the 
final denoised map). The minimization is finally achieved in 
the wavelet domain by choosing a threshold value that min- 
imizes the risk at each wavelet level (see for instance Ogden 
1997). 

Results of the errors in the map reconstruction are 
shown in Table 1. The map error is defined as: 
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(21) 



Performing 20 simulations (proved to be enough as results 
reached stable values) at each S/N level we have also cal- 
culated the lcr error. The error improvement achieved with 
the denoising technique applied goes from factors of 3 to 5 
for S/N = 3 to S/N = 0.7. 

It is also interesting to see how well the denoising 
method performs to reconstruct the temperature power 
spectrum. Mean values and lcr errors of the relative er- 
rors, \ACt/Ce\, axe shown in Figure 3 for three S/N ratios 
and the case of non-uniform noise considered in this work. 
The C/s are reconstructed from the denoised maps with 
\AC e /Ce\ < 10% up to I ~ 1000 in cases S/N > 1. This er- 
ror can only be achieved up to an I ~ 700 in the S/N — 0.7 
case. Higher order multipoles (i ^ 1500) are reconstructed 



Table 1. Reconstruction errors vs S/N. 



S/N 


% 


map error ±lcr 


0.7 




26.3±0.4 


1.0 




20.7±0.4 


2.0 




13.3±0.2 


2.0 (n.u.) 




14.3±0.3 


3.0 




10.3±0.2 



with \ACt/Ce\ ^ 20%. Absolute relative errors and recon- 
structed Ci 's for a given map are presented in Figures 4 and 
5 respectively. 

In order to check the performance of the SURE thresh- 
olding technique, knowing the original maps we can find 



© 1998 RAS, MNRAS 000, §-?? 



Table 2. Reconstruction errors vs threshold, S/N=l. 



Threshold 


% map error 


hard 


23.5 


1.5 Un 


21.7 


1.0 cr n 


20.7 


0.7 cr n 


20.5 


0.6 a„ 


20.6 


0.5 CT n 


20.7 


0.4 ff„ 


21.0 


0.3 cr n 


21.3 


signal+noise 


100.0 



the optimal threshold to get a reconstructed map with a 
minimum error (as defined above). In S/N = 1 maps the 
optimal threshold is found to be 0.6 — 0.7a n . Thresholds 
between 0.3 — la„ do not make substantial changes in the 
reconstructed map (see Table 2). The hard case included in 
that table stands for a case where all coefficients below a 
signal-to-noise dispersion ration < 1.5 are removed, leaving 
the others untouched. For comparison, the error obtained 
comparing the signal plus noise map with the original signal 
map is also presented in Table 2. We can see that the error 
reconstruction achieved with the SURE technique equals the 
one obtained with the optimal threshold. 

A comparison of wavelet techniques with Wiener fil- 
ter (see for instance Press et al. 1994) has also been per- 
formed. In relation to map reconstruction the error affecting 
the Wiener reconstructed maps is comparable to the error 
for the wavelet reconstructed maps, in all cases. However, 
in order to apply Wiener filter previous knowledge of sig- 
nal power spectrum is requiered. Reconstructed and resid- 
ual maps using both, wavelets and Wiener filter, are shown 
in Figure 6. Regarding the Ces, performance of Wiener fil- 
ter is clearly worse than Wavelets for £ = 1000 — 1500, as 
can be seen in Figure 4. For example, for a S/N = 1 the 
Ces are recovered using Wiener filter with an error between 
20% and 70% for is between 1000 and 1500 being this error 
smaller by a factor of 2-4 for Wavelet reconstruction. The 
error is also clearly larger for Wiener reconstruction than 
for Wavelet reconstruction, up to £ ~ 2000 in cases with 
S/N > 1. 

We have checked for non-Gaussian features possibly in- 
troduced by the non-linearity of the soft thresholding used 
in the wavelet methods applied for denoising. Distributions 
of Skewness and Kurtosis have been obtained for the original 
signal maps as well as for the denoised ones. No significant 
differences can be appreciated between both distributions. 
However this method could not be good enough to detect 
non-Gaussian features. As recently claimed by Hobson et 
al. (1998), the analysis of the distribution of wavelet coeffi- 
cients is one of the most efficient methods to detect them. 
We have performed a similar analysis using the Daubechies 
4 multiresolution wavelet coefficients. These coefficients are 
Gaussian distributed in the case of a temperature Gaus- 
sian random field. The application of soft thresholding to 
the wavelet coefficients at a certain level clearly changes the 
Gaussian distribution by removing all coefficients whose ab- 
solute values are below the imposed threshold and shifting 



Figure 6. 12°. 8 X 12°. 8 maps of the cosmological signal (top 
left), signal plus noise with S/N = 1 (top right), denoised map 
using a soft thresholding as explained in the text (middle left) 
and residual map obtained from the CMB signal map minus the 
denoised one (middle right). For comparison a denoised map using 
Wiener filter is presented in the bottom left panel together with 
the residuals in the bottom right panel. 



the remaining ones by that threshold. As an example, in the 
previous case S/N — 1 the kurtosis of the diagonal level 
3 distribution changes from 3.3 ± 0.1 to 34 ± 10 ! (notice 
that the change strongly depends on the threshold imposed) . 
This result is not surprising as any non-linear method used 
for denoising or foreground separation will introduce non- 
Gaussinity at different levels in the reconstructed map. For- 
tunately there are two ways of overcoming the question of 
determining the Gaussianity of the CMB signal. One way 
would be to check the Gaussian character of the data before 
applying denoising to maps affected by Gaussian noise. We 
have checked this by looking at the multiresolution wavelet 
coefficient distributions in the case of S/N = 1. The addition 
of white noise didn't change the mean value and error bar 
of the kurtosis. The second way would be applying a linear 
denoising method. We have used a simple one consisting in 
removing all detail coefficients at levels with signal-to-noise 
dispersion ratio < 1.5 (notice that 1.5 corresponds the upper 
value of the threshold interval where soft thresholding was 
applied) . This method is equivalent to applying hard thresh- 
olding with a threshold above all the coefficients. The errors 
of the reconstructed map and its corresponding C;s increase 
slightly compared to the SURE thresholding method (see 
table 2and top- left panel of figure 4). The same hard thresh- 
olding linear method will give even better results using 2-D 
Wavelets with two scales of dilation (Sanz et al. 1999) in- 
stead of the one-scale multiresolution techniques, since the 
former works with many more resolution levels being there- 
fore more selective in removing the coefficients. 



5 DISCUSSION AND CONCLUSIONS 

A wavelet multiresolution technique has been presented and 
used to analyse and denoise CMB maps. This method has 
been proved to be one of the best to reconstruct observed 
CMB maps as well as power spectra by removing a signifi- 
cant percentage of the noise. The analysis has been carried 
out assuming a uniform Gaussian noise as would be expected 
in a small sky patch, e.g. 12°. 8 x 12°. 8, observed by satel- 
lite scans. Analysis of whole sky CMB maps using wavelets 
will be performed in a future work. Since these data will be 
affected by non-uniform noise, the use of wavelet techniques 
to localize map features will be highly suitable. 

A semi-analytical calculation of the variance of the 
wavelet coefficients has been presented. The behaviour of 
the variance of the detail coefficients is given for a standard 
CDM model in the case of Haar and Mexican Hat bases. 
The acoustic peaks can be noticed in the wavelet coefficient 
variance represented in Figure 1. Moreover, these peaks are 
better defined for the Mexican Hat wavelet system since 
these wavelets are more localized than the Haar ones. 

Denoising of CMB maps has been carried out by using 
a signal-independent prescription, the SURE thresholding 
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method. The results are model independent depending only 
on the observed data. However, a good knowledge of the 
noise affecting the observed CMB maps is required. For a 
typical case of S/N ~ 1 the high order detail coefficients are 
dominated by the signal, whereas the lowest ones are noise 
dominated. This behaviour is due to the expected depen- 
dence of the temperature power spectrum, Ce oc l~ 2 . The 
applied wavelet method is able to reconstruct maps with an 
error improvement factor between 3 and 5 and the CMB 
power spectrum of the denoised maps carries relative errors 
below 20% up to / ~ 1500 for S/N > 1. We have also checked 
that SURE thresholding methods are providing thresholds 
in agreement with the optimal ones. 

For comparison Wiener filter has also been applied to 
the simulations considered in this paper. This method recon- 
structs CMB maps after denoising with errors comparable to 
the Wavelet method we propose, as shown in Figure 6. How- 
ever, the C7s of the denoised maps obtained applying Wiener 
filter have relative errors larger than a factor of 2 than the 
relative errors of the Ces obtained from the wavelet recon- 
structed maps in the range £ = 1000 — 1700. In addition we 
have applied a Maximum Entropy Method (MEM) to the 
maps used in this work, with the definition of entropy given 
by Hobson & Lasenby (1998). This method provides recon- 
struction errors at the same level as multiresolution wavelet 
methods. However, the later are easier (not requiring itera- 
tive processes) and faster (0(N)) to apply than MEM. 

A possible handicap of denoising methods based on soft 
thresholding of wavelet coefficients as well as other non- 
linear methods are the non-Gaussian features introduced in 
the reconstructed map. However one can still detect the pos- 
sible intrinsic non-Gaussianity of the CMB signal by study- 
ing it in the signal plus noise map using the wavelet coef- 
ficient distribution. Moreover a valid reconstruction can be 
obtained by applying a "hard" thresholding linear method 
as discussed in the text. 

In a different work, we are studying the case of using 
a wavelet method based on two scales of dilation (Sanz et 
al. 1999). Though this method has the advantage of keep- 
ing information on two different scales, for the purpose of 
denoising both methods give comparable results. The lin- 
ear hard thresholding method is expected to perform better 
for 2-D wavelets than for multiresolution ones as the former 
works with many more resolution levels. 

Summarizing, the main advantages of the wavelet 
method are: to provide local information of the contribu- 
tion from different scales, to be computationally very fast 
0(N), absence of tunning parameters and the most impor- 
tant, the good performance on denoising CMB maps. The 
best reconstruction is achieved using soft thresholding tech- 
niques. Concerning the Gaussianity of the signal one can 
apply the suggested linear method for denoising. Moreover, 
the soft thresholding technique will provide a good reference 
map and power spectrum for the signal, that can be used 
to check the quality of other reconstructions based on linear 
methods. Wavelets are also expected to be a very valuable 
tool to analyse future CMB maps as those that will be pro- 
vided by future missions like MAP and Planck. 
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